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ABSTRACT 

Two small satellites of Pluto, S/2005 PI (hereafter PI) and S/2005 P2 (hereafter 
P2), have recently been discovered outside the orbit of Charon, and their orbits are 
nearly circular and nearly coplanar with that of Charon. Because the mass ratio of 
Charon-Pluto is ~ 0.1, the orbits of P2 and PI are significantly non-Keplerian even 
if P2 and PI have negligible masses. We present an analytic theory, with P2 and 
PI treated as test particles, which shows that the motion can be represented by the 
superposition of the circular motion of a guiding center, the forced oscillations due to the 
non-axisymmetric components of the potential rotating at the mean motion of Pluto- 
Charon, the epicyclic motion, and the vertical motion. The analytic theory shows that 
the azimuthal periods of P2 and PI are shorter than the Keplerian orbital periods, and 
this deviation from Kepler's third law is already detected in the unperturbed Keplerian 
fit of Buie and coworkers. In this analytic theory, the periapse and ascending node of 
each of the small satellites precess at nearly equal rates in opposite directions. 

From direct numerical orbit integrations, we show the increasing influence of the 
proximity of P2 and PI to the 3:2 mean-motion commensurability on their orbital 
motion as their masses increase within the ranges allowed by the albedo uncertainties. 
If the geometric albedos of P2 and PI are high and of order of that of Charon, the 
masses of P2 and PI are sufficiently low that their orbits are well described by the 
analytic theory. The variation in the orbital radius of P2 due to the forced oscillations 
is comparable in magnitude to that due to the best-fit Keplerian eccentricity, and there 
is at present no evidence that P2 has any significant epicyclic eccentricity. However, 
the orbit of PI has a significant epicyclic eccentricity, and the prograde precession of 
its longitude of periapse with a period of 5300 days should be easily detectable. If the 
albedos of P2 and PI are as low as that of comets, the large inferred masses induce 
significant short-term variations in the epicyclic eccentricities and/or periapse longitudes 
on the 400~500-day timescales due to the proximity to the 3:2 commensurability. In 
fact, for the maximum inferred masses, P2 and PI may be in the 3:2 mean-motion 
resonance, with the resonance variable involving the periapse longitude of PI librating. 
Observations that sample the orbits of P2 and PI well on the 400-500-day timescales 
should provide strong constraints on the masses of P2 and PI in the near future. 
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1. INTRODUCTION 

Weaver et al. (2006) have recently discovered two new satellites of Pluto, S/2005 PI (hereafter 
PI) and S/2005 P2 (hereafter P2), from images taken with the Hubble Space Telescope (HST). These 
are the first new satellites of Pluto since the discovery of Charon in 1978 (Christy and Harrington 
1978). The new satellites are much fainter and hence much smaller than Charon (whose diameter 
~ 1200km; e.g., Sicardy et al. 2006), with the diameter of PI 60-170 km depending on the 
geometric albedo, and P2 about 20% smaller than PI. Since the discovery data consist of only 
two epochs separated by three days in May 2005, Weaver et al. (2006) were unable to determine 
the orbits of P2 and PI, but the data are consistent with the orbits of P2 and PI being nearly 
circular and nearly coplanar with that of Pluto-Charon, with orbital periods of ^ 25 days for P2 
and ~ 38 days for PI. Steffi et al. (2006) have used the same HST observations to place constraints 
on the existence of any additional satellites of Pluto. 

Prior to the discovery of P2 and PI, Buie et al. (2006, hereafter BGYYS) have obtained HST 
images of the Pluto system in a series of 12 visits from June 2002 to June 2003 for the purpose of 
producing an albedo map of Pluto. BGYYS were able to detect Charon in individual frames and 
P2 and PI by stacking the images taken during each visit. BGYYS fit this data, together with 
the data for Charon from Tholen and Buie (1997) and the data for P2 and PI from Weaver et 
al. (2006), by assuming that all three satellites are on unperturbed Keplerian orbits. The best-fit 
Keplerian orbital parameters found by BGYYS and their 1-a errors are reproduced in Table 1. The 
best-fit orbit of Charon relative to Pluto is consistent with zero eccentricity, and BGYYS pointed 
out that the nonzero eccentricity reported by Tholen and Buie (1997) is probably due to the use 
of an imprecise center of body of Pluto in the earlier paper. The parameters in Table 1 show that 
the orbits of P2 and PI are indeed nearly circular and nearly coplanar with that of Pluto-Charon. 
Because P2 and PI are much smaller than Pluto and Charon, they orbit about a point that is very 
close to the center of mass of Pluto and Charon. Thus one of the fitting parameters is the mass 
ratio of Charon-Pluto, and BGYYS found mc/nip = 0.1165 ± 0.0055. 

Two results from the orbit-fitting suggest that unperturbed Keplerian orbits are not good 
assumptions for the orbits of P2 and PI. The orbital period P and semimajor axis a of each orbit 
are independent parameters in the fits by BGYYS. If we assume that Kepler's third law is valid for 
all of the orbits, each orbit yields an independent measurement of the total mass of Pluto-Charon. 
BGYYS found rup + rric = 1.4570 ± 0.0009 x 10^^ kg from Charon's orbit, 1.480 ± 0.011 x 10^^ kg 
from P2's orbit, and 1.4765 ± 0.006 x 10^^ kg from Pi's orbit. ^ The values of nip + rric from P2 and 
PI disagree with that from Charon at the 2.1-cj and 3.2-0" levels, respectively. As we shall see, the 
discrepancies are in fact due to the deviation from Kepler's third law for P2 and PI. Because of 
the rather large mass ratio of Charon-Pluto, the deviations of the gravitational potential from that 
of a point mass and, consequently, of the orbits of P2 and PI from Keplerian orbits are nontrivial. 



We note that the masses reported by BGYYS are consistent with using G = 6.67 X 10"" kg-^ s"^ instead of 
G = 6.672 X 10~" kg~^ s~^ in converting from G{mp + nic) = (27r/P)^a^ to rUp + rric. 
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even if P2 and PI can be treated as test particles. 

The second result from the orbit-fitting that suggests non-Keplerian orbits for P2 and PI 
is the confirmation of the result of Weaver et al. (2006) that the orbital periods of Charon, P2, 
and PI are nearly in the ratio 1:4:6. This means that the orbits of P2 and PI could be strongly 
affected by resonant or near-resonant interactions. As we shall see, the strongest effects come from 
the proximity of P2 and PI to the 3:2 mean-motion commensurability, which is the lowest order 
commensurability among the satellites, even though P2 and PI are much smaller than Charon. 

In Section 2 we present an analytic theory for the orbits of P2 and PI that is valid in the 
limit that the satellites have negligible masses and can be treated as test particles. It shows that 
the motion can be represented by the superposition of the circular motion of a guiding center, 
the forced oscillations due to the non-axisymmetric components of the potential rotating at the 
mean motion of Pluto-Charon, the epicyclic motion, and the vertical motion. It also gives analytic 
results for the deviation from Kepler's third law and the periapse and nodal precession rates. 
In Section 3 we present direct numerical orbit integrations with different assumed masses for P2 
and PI within the ranges allowed by the uncertainties in the albedos. The numerical results are 
compared to the analytic theory in Section 2, and the increasing importance of the proximity to 
the 3:2 commensurability with increasing masses is examined. In fact, for the maximum masses 
corresponding to the lowest expected albedos, P2 and PI may be in the 3:2 mean-motion resonance, 
with the resonance variable involving the periapse longitude of PI librating. In Section 4 we 
summarize our results and discuss the prospects for detecting non-Keplerian behaviors and putting 
constraints on the masses of P2 and PI with existing and future observations. 



2. ANALYTIC THEORY 

In this section we develop an analytic theory for the orbits of the satellites P2 and PI that is 
valid in the limit that the satellites have negligible masses and can be treated as test particles. The 
latest orbital fit show that the orbit of Charon relative to Pluto is consistent with zero eccentricity 
(see Table 1). Thus we assume that the orbit of Charon relative to Pluto is Keplerian and circular, 
with semimajor axis a.pc and mean motion (or circular frequency) Upc = [G{mp + rric) / a^(\^^'^ 
where nip and rUc are the masses of Pluto and Charon, respectively. In a cylindrical coordinate 
system with the origin at the center of mass of the Pluto-Charon system and the z = plane 
being the orbital plane of Pluto-Charon, the positions of Charon and Pluto are Vc = {ac,(t>ci^)i 
and Tp = {ap,(f)c + 7r,0), respectively, where ap = apcmd {jUp + mc), ac = apcvrip/ (rup + rric), 
4>c{t) = npct + ippc, and ippc is a constant. 
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2.1. Potential 

Since the orbital radii of the sateUites are much smaUer than the Hih radius (~ 8.0 x 10^ km) 
of the Pluto-Charon system, the perturbation from the Sun can be ignored and the gravitational 
potential at r = {R, cf), z) is 



Grrip Gnic 
$(r) = —r - 



(1) 



The orbits of the satellites are nearly coplanar with that of Pluto-Charon, and we expand l/jr" — rd 
in powers of z: 



where 



p= [R^ + al- 2Rac cos ((?!) - (t)c)] . 



(3) 



By expressing the inverse powers of p as cosine series using the Laplace coefficients (see Eq. [6.62] 
of Murray and Dermott 1999), we obtain 
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where Sko is the Kronecker delta, etc = cic/R, and 6j(ac) are the Laplace coefficients. With a similar 
expression for l/|r — rp|, the potential can be written as 
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The axisymmetric fc = components of the potential are identical to those due to two rings — one 
of mass nip and radius ap and another of mass nic and radius Cc — at the z = plane. 

Since apc/R ~ 0.40 and 0.30 for P2 and PI, respectively, it is instructive to examine the 
expansion of <I>ofc in powers of apc/R for the lowest values of k: 



00 



1 + 



1 



4(1 -|- nic/nip)'^ 



VmJ V R ) 



(7) 



+ 



9{1 — mc/nip + nic /nip) f 



nir 



64(1 + nic/nip)^ V^p/ ^ 



apc\ 4 



— ^ + 



G{mp + nic 
R 



01 



3{l — nic/nip) f nic 
S{1 + nic/nipY \nT'pJ ^ R 



+ •■ 



G{nip + nic 
R 



(8) 



-5- 



02 



m„ J \ R / 



4(1 + mc/mp) y/<<,p 

5(1 - rric/mp + ml/m^) / m, 



(9) 



16(1 + rric/mpY 



m 



+ 



G{mp + rric 
R 



03 



5(1 - irijni^) 
8(1 + rric/mp)^ 

and also the expansion of $20: 
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Prom the series definition of the Laplace coefficients (Eq. [6.68] of Murray and Dermott 1999), one 
expects the lowest order terms of to be of the order of {apc/R)^. However, the terms of the 
order of Upc/R in $01 cancel, and the lowest order nonzero terms in $01 is of the order of {upc/R)^. 
Eqs. (7) and (11) show that the axisymmetric components of the potential deviate from that of 
a point mass of mass nip + rric at the origin by terms of the order of {apc/R}^ and higher, while 
Eqs. (8)-(10) show that the non-axisymmetric components of the potential are weak and of the 
order of (apc/R)^ and higher. It should be noted that the deviations from a point mass potential 
are multiplied by an additional small quantity rric/mp. 



2.2. Equations of Motion and Solution 

The equations of motion in cylindrical coordinates are 

.. _ 5$ 

dz 

The gravitational potential (Eq. [5]) is weakly non-axisymmetric and rotates at a constant pattern 
speed Upc, and the satellites P2 and PI are on nearly circular orbits that are nearly coplanar with 
that of Pluto-Charon. An approximate solution to the equations of motion that is valid for such 
orbits can be obtained by representing the orbit as small deviations from the circular motion of a 
guiding center in the z = plane: 

R = Rq + Ri{t), 

(t> = Mt) + Mt), (13) 

Z = Zl{t), 
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where the constant Rq is the radius of the guiding center, \Ri/Ro\ <C 1, |^i| <C 1, and |zi/i?o| <S 1- 
The approach is well known from the theory of orbits in weakly barred galaxies (sec, e.g., Binney 
and Tremaine 1987) and planetary ring dynamics (see, e.g., Goldreich and Tremaine 1982). 



Substituting Eqs. (5) and (13) into Eq. (12), the only nontrivial equation at the zeroth order 
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which describes the circular motion of the guiding center. The solution of Eq. (14) is 

4>o{t) = not + (po, 
where ipo is a constant and the mean motion no is given by 
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In the above equation nx = [G{mp + nic) / R^^^"^ is the Keplerian mean motion at Rq, D = d/da, 
and Oc and Qp are evaluated at R = Rq. 

At the first order, the equation for is 

7 2no • ^fc$ofe(-Ro) . n AA\ no^ 
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An integration of Eq. (18) yields an expression for (pi, which can then be substituted into the first 
order equation for Ri to yield 
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where the epicyclic frequency kq is given by 
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In Eqs. (20) and (21), n = {R~^d^oo/dRy^^ is the mean motion at R, and the quantities in the 
square brackets are evaluated at R = Rq. Eq. (20) is the equation of motion of a simple harmonic 
oscillator of natural frequency kq that is driven at frequencies k\no — n^d- It is straightforward to 
solve Eqs. (18) and (20) to obtain 
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Finally, for the motion in z, the first order equation for zi is 



Zl + VqZi = 0, 



and its solution gives 



Z = Zl = RqI cos(i'of + C) 
where i and C are constants and the vertical frequency uq is defined by 
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The motion in R and (j) is described by the superposition of the circular motion of the guiding 
center at Rq at frequency uq, the epicychc motion represented by "eccentricity" e at frequency kq, 
and the forced oscillations of fractional radial amplitudes Ck at frequencies k\nQ — npc\. The motion 
in z decouples from that in R and (p and has only free oscillations at the vertical frequency i/q- If 
there is no epicyclic or vertical motion [e = i = 0), the orbit is a closed orbit in the frame rotating 
at the pattern speed Upc, with the variation of the radius with the azimuthal angle being the sum 
of terms with fractional amplitude Ck and period 2'!r/k. If e S> Ylk the orbit is approximately 
a precessing Keplerian ellipse with eccentricity e, mean motion no, and periapse precession rate 

w = no — kq. (33) 

If in addition i ^ 0, the orbit has an inclination i with respect to the Pluto-Charon orbital plane 
and a nodal precession rate 

J7 = no — t'o- (34) 

2.3. Deviation from Kepler's Third Law 

Eq. (17) shows that the mean motion no and hence the azimuthal period Pq = 27r/no deviate 
from those of a Keplerian orbit. In Fig. 1 the solid line shows Pr/Pq — 1 = no/riK — 1 from 
Eq. (17) as a function of the guiding center radius Ro for the best-fit mass ratio of Charon-Pluto, 
rric/rnp = 0.1165, while the dotted lines show the same quantity for rric/rnp that arc 1 a (±0.0055) 
from the best-fit value. (The uncertainty in mc/mp dominates that in ape in the evaluation of Eq. 
[17].) Fig. 1 shows that the azimuthal period is shorter than the Keplerian orbital period. 

The orbital parameters of P2 and PI in Table 1 were obtained by BGYYS from unperturbed 
Keplerian fits. As we have shown in this section, the orbits of P2 and PI are non-Keplerian even if 
they can be treated as test particles. Nevertheless, the orbital period P and semimajor axis a (which 
are independent parameters in the fits by BGYYS) in Table 1 should closely resemble the azimuthal 
period Pq and guiding center radius Rq. With Pk = 2n[a^ /G{mp -|-TOc)]^/^ = Ppc{o./ CLpcf^"^ , where 
Ppc and flpc are the orbital period and semimajor axis of Pluto-Charon from Table 1, we find 
Pk/P - I = 0.0079 ± 0.0038 and 0.0067 ± 0.0021 for P2 and PI, respectively The values of 
Pk/P — 1 and a for P2 and PI are shown in Fig. 1 with their 1 a error bars. The orbital periods 
are clearly shorter than the Keplerian values (by 2.1a and 3.2 cr for P2 and PI, respectively). On 
the other hand, Pr/P — 1 for P2 is in excellent agreement (within 0.4(7) with the analytic result, 
and that for PI is in reasonable agreement (within 1.6 a") with the analytic result. The remaining 
discrepancy for PI may be simply statistical, but it could also be due to the assumption in the 
fitting that the orbit is an unperturbed Keplerian orbit (see below for more details on the expected 
non-Keplerian behaviors). 
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2.4. Periapse and Nodal Precessions 

With a from Table 1 as i?o for P2 and PI and rric/mp = 0.1165, we evaluate Pk = Ppc{C'/apc)^^'^ , 
no/riK (Eq. [17]), kq/uk (Eq. [22]), and vqItik (Eq. [32]) for the analytic theory, and they are listed 
in Table 2. The precession of the periapse is prograde with period li^ l\w\ = 27r/|no — ko\ = 1740 
and 5280 days for P2 and PI, respectively. The nodal precession has a similar period (27r/|0| = 
27r/|no — i^ol = 1770 and 5330 days for P2 and PI, respectively) but it is retrograde. 

The periapse and nodal precessions at nearly equal rates in opposite directions and the faster- 

than-Kcplerian mean motion are similar to the behaviors of orbits around an oblate planet (see, 
e.g.. Section 6.11 of Murray and Dermott 1999). This can be understood from the fact that the 
{ttpc/Kf' terms in the axisymmetric components of the potential, <I>oo and $20 in Eqs. (7) and (11), 
are identical to the J2 terms of an oblate planet with J2 = mpmc/[2(mp + mc)^]. 

3. NUMERICAL ORBIT INTEGRATIONS 

3.1. Initial Conditions and Numerical Methods 

For the remainder of this paper we use Jacobi coordinates where the position of Charon is 
relative to Pluto, the position of the inner satellite P2 is relative to the center of mass of Pluto- 
Charon, and the position of the outer satellite PI is relative to the center of mass of Pluto-Charon- 
P2. Jacobi coordinates are the natural generalization of the coordinates used in Section 2 (where 
the position of PI is relative to the center of mass of Pluto-Charon) when P2 and PI are not test 
particles, and they reduce to the coordinates used in Section 2 in the test-particle limit. 

Prom Ppc and a^c in Table 1, we adopt G{mp + rric) = {i-K/Ppcfa^^ = 9.71791 x 10^^ m^s~^ 
(or rup + mc = 1.4565 x 10^^ kg for G = 6.672 x 10~^^ m^ kg^^ s~^). For the mass ratio rric/mp, we 
use the best-fit value 0.1165 from BGYYS. We generate the initial position and velocity of Charon 
relative to Pluto by using the orbital parameters in Table 1 at epoch JD 2452600.5 as the osculating 
Keplerian orbital parameters. 

The orbits of P2 and PI arc sufficiently non-Keplerian even in the test-particle limit that, if we 
had generated their initial conditions by assuming that the orbital parameters in Table 1 are the 
osculating Keplerian parameters at epoch JD 2452600.5, the numerically integrated orbits would 
have properties significantly different from those of the best-fit Keplerian orbits in Table 1, with 
the mean orbital radii smaller than the semimajor axes in Table 1 and the variations in the orbital 
radii larger than those for the eccentricities in Tabic 1 . Without refitting the data, wc aim to adopt 
a set of initial conditions so that the resulting orbits would closely resemble the best-fit Keplerian 
orbits in Table 1. This is accomplished by using the mean longitudes A of P2 and PI in Table 1 
as the initial values of ^ and then using Eqs. (23), (25), and (26) to set the initial values of R, R, 
and <j). As in Sections 2.3-2.4, we adopt the semimajor axes a from the Keplerian fit as the guiding 
center radii Rq. We include the forced oscillation terms up to A; = 4, and the coefficients Ck and 
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Dk are listed in Table 2, along with Pk-, riQ/riK, no/riK-, and vo/tik- For the initial values of 
in the forced oscillation terms, we can ignore the difference between <j) and (f)Q (with the latter for 
the guiding center) and use ^ — ^c- 

We can see from Table 2 that X^^, Ck ~ 0.0029 for P2. Thus the fractional orbital radius 
variation due to the forced oscillation terms alone is comparable to that due to the best-fit Kcplerian 
eccentricity (0.0023), which is itself consistent with zero (±0.0021). Therefore, we set the initial 
epicyclic eccentricity e = for P2. For PI, ~ 0.0004, which is much smaller than the best-fit 

Keplerian eccentricity of 0.0052(ib0.0011), and we adopt an initial epicyclic eccentricity e = 0.0052. 
In the epicyclic approximation for a Keplerian orbit, the phase kqI + ^ is the mean anomaly. So it 
is reasonable to adopt for PI i/j = X — w, where A and w are the mean longitude and longitude of 
periapse from Table 1. 

BGYYS have computed the 1-a contours of the orbit poles on the J2000 sky plane for their 
best-fit Keplerian orbits for Charon, P2, and PI using the Monte Carlo technique, and the contours 
are shown in their Fig. 2. The 1-a contour for Charon is significantly smaller than those for P2 
and PI, with the latter two having mean radii of about 0?34 and 0?16, respectively. The orbit 
pole of Charon is 0?10 from that of P2 for the best-fit orbits, which is well within the 1-a contour 
for P2, and it is 0?25 from that of PI for the best-fit orbits, which is about 50% further than 
the 1-a contour for PI and only marginally significant. Since there is no significant detection of 
any mutual orbital inclinations, we assume that the orbits of P2 and PI are coplanar with that 
of Pluto-Charon. This means that the precession of nodes is not examined by our numerical orbit 
integrations. 

We perform 5 sets of integrations with different assumed masses for P2 and PI. From the 
photometry of P2 and PI, Weaver et al. (2006) have estimated that the diameters of P2 and PI 
are 46 it 4 km and 61 it 4 km, respectively, if the geometric albedos are Charon-like and = 0.35. On 
the other hand, if the albedos are comet-like and = 0.04, the diameters are 137 ± 11 km for P2 and 
167 lb 10km for PI. If we assume that the mean density is 2gcm~^ (i-e., similar to that of Pluto), 
771-2 = 1-02 X 10^'' kg and rrii = 2.38 x 10^'' kg in the high albedo case, and m2 = 2.69 x 10^^ kg 
and mi = 4.88 x 10^^ kg in the low albedo case. In addition to two integrations with the high- and 
low-albedo masses, we perform an integration with masses 10~^ times those of the high albedo case 
(so that P2 and PI are test particles), an integration with masses twice those of the high albedo 
case, and an integration with masses half those of the low albedo case. 

The direct numerical orbit integrations are performed using a modified version of the Wisdom 
and Holman (1991) symplectic integrator contained in the SWIFT^ software package. The Wisdom- 
Holman integrator is based on dividing the Hamiltonian of the gravitational iV-body problem into 
a part that describes the Keplerian motions of the satellites around the planet (or, in the case of 
a planetary system, the planets around the star) and a part that describes the perturbations to 



^See http://www.boulder.swri.edu/~hal/swift.html. 
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the Keplerian motions. The division used by Wisdom and Holman (1991) assumes that all of the 

satellite masses are much smaller than the planet mass, but in the case of Pluto, rric/mp = 0.1165. 
We have described in Lee and Peale (2003) a modified Wisdoni-Holnian integrator using a slightly 
different division of the Hamiltonian into the Keplerian and perturbation parts. The modified 
integrator is designed for hierarchical systems, where the masses of the satellites can be comparable 
to that of the planet but the orbit of each satellite is much larger than that of the satellite just inside, 
and it was used by Lee and Peale (2003) to study hierarchical extrasolar planetary systems. This 
modified integrator can also handle the Pluto system, where {mc/'mp){apc/ R)^ ^ 0.019, without 
an excessively small timcstep. The integrations are performed with a timestep of 10^ s (or about 
55 steps per Charon's orbit). 



The numerical orbit integrations with test-particle masses, high-albedo masses, and twice the 
high-albedo masses for P2 and PI show that the orbits of P2 and PI are well described by the 
analytic theory in Section 2 if the masses of P2 and PI are of the order of the high albedo ones. In 
Fig. 2 we plot the variations in the orbital radii R2 and i?i of P2 and PI, respectively, for 800 days 
in the high albedo case. The dashed lines indicate the semimajor axes, a, and the maximum and 
minimum radii, a(l ± e), for the best-fit Keplerian orbits in Table 1. The orbital radius R2 of P2 
clearly shows the forced oscillations, which are dominated by the Ci and C2 terms (see Table 2) 
with periods 27r/|no — n^d ~ 8.6 days and 7r/|no — npc\ ~ 4.3 days, respectively. As expected, the 
variation in R2 due to the forced oscillations is comparable in magnitude to that due to the best-fit 
Keplerian eccentricity. Although the initial epicyclic eccentricity of P2 is set to zero, the fact that 
the forced oscillation terms with k > A are not included in setting the initial conditions results in 
a small epicyclic motion with a period 27r/Ko ~ 25.2 days, which is also visible in the plot of R2 
in Fig. 2. The variation in R\ is dominated by the epicyclic motion with e « 0.0052 and period 
27r/Ko ~ 38.6 days. For PI, the amplitudes of the forced oscillations are significantly smaller than 
for P2, as we expect from the values of Cfe in Table 2, and the dominant forced oscillation periods 
are 27r/|no — ripd ~ 7.7 days and 7r/|no — ripd ~ 3.8 days. The fc > 4 terms are sufficiently small for 
PI that their neglect in the initial conditions do not result in any significant additional epicyclic 
motion. 

In order to study in more detail the epicyclic component of the motion, we need to eliminate 
the high frequency forced oscillations. We find that the forced oscillations are sufficiently close 
to those predicted by the analytic theory in all of our numerical integrations that they can be 
effectively eliminated by defining a transformed orbital radius: 



k 

with Rq and C}~ from Table 2 and (j) and from the numerical integrations themselves (compare 
to Eq. [23]). In Fig. 3 we plot the variations in R'2 and R'^ in the high albedo case, and it is clear 



3.2. 



Results 




(35) 
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from a comparison with Fig. 2 that most of the forced oscillations are eliminated. Fig. 3 shows 
that there are small periodic variations in the maximum and minimum values of R'2 and i?^, or 
equivalently in the epicyclic eccentricities 62 and ei, in this run with the high-albedo masses. The 
amplitudes of the eccentricity variations are larger in the run with twice the high-albedo masses, 
while the eccentricities are nearly constant in the run with test-particle masses, which indicate that 
the variations are due to interactions between P2 and PI (see below for more details). 

The azimuthal period Pq can be determined for the numerical integrations from the cumulative 
increase in <j). We find Pq = 24.913 and 38.335 days for P2 and PI, respectively, for the run with 
high-albedo masses, and they are identical to Pq found for the run with test-particle masses.^ If 
we use Pk from Table 2, no/riK = Pk/Pq = 1-0056 and 1.0033 for P2 and PI, respectively. The 
latter is in excellent agreement with the analytic value in Table 2, but the former is slightly smaller 
than the analytic value. The small discrepancy for P2 is due to the fact that the numerically 
integrated orbit has a non-zero epicyclic eccentricity, with P2 near the periapse at t = (see Fig. 
3) , which means that the guiding center radius of the numerically integrated orbit [Rq = 48698 km 
from the average of the maximum and minimum values of -R2) is slightly larger than the value we 
were aiming for [Rq = 48675 km). (For PI, the guiding center radius of the numerically integrated 
orbit is identical to the value we were aiming for.) If we use Pk = 25.0695 days for the numerically 
determined Rq instead of P^: = 25.0518 days from Table 2 for P2, we find riQ/riK = Pk/Pq = 1.0063, 
which is in excellent agreement with the analytic value in Table 2 (note that the ratio no/nx is 
much less sensitive to a small change in Rq than Pk and Pq separately). 

The easiest way to determine the longitude of periapse w of the epicyclic motion as a function 
of time is to monitor the transformed orbital radius R' (Eq. [35]) during the numerical orbit 
integration. When R' changes from decreasing in a previous timestep to increasing in the current 

timestep, the satellite has passed the periapse, and we can use the values of R' and (p at the end 
of the current step and two previous steps to find the time of periapse passage and ru {= (p ai 
the time of periapse passage) by interpolation. Similarly, the time at which w = (p + 180° (i.e., 
apoapse passage) can be determined. Fig. 4 shows the evolution of tU2 and wi of P2 and PI, 
respectively, for 10^ days in the high albedo case. There are some spurious points in the plot of 
W2, because the forced oscillations arc not perfectly eliminated in the transformed orbital radius 
and there are occasional false minima and maxima in i?2 when the epicyclic eccentricity 62 is very 
small (see Fig. 3). In this high albedo case, the long-term evolution of TU2 and zui are prograde 
precessions with periods of 2000 and 5300 days, respectively. The latter is in excellent agreement 
with the analytic result in Section 2.4, but the former is about 15% longer. The 15% discrepancy 



^It should be recalled that we adopt the best-fit semi-major axes a from the Keplerian fit as the guiding center 

radii _Ro and that the 1-a error in a is ±121 km for P2 and ±88 km for PI (see Table 1). If we vary the guiding 
center radii by ±1 a in the initial conditions of our numerical integrations, we would find Po = 24.913 ± 0.093 and 
38.335 ± 0.078 days for P2 and PI, respectively. Thus, similar to the analytic results in Section 2.3, the azimuthal 
periods from the numerical integrations agree with the best-fit orbital periods in Table 1 from the Keplerian fit to 
within 0.6 a for P2 and 1.6 a for PI. 
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for P2 is too large to be explained by the slightly larger guiding center radius mentioned in the 
previous paragraph, and it is not due to interactions between P2 and PI, as the tcst-particic run 
shows the same precession period. We suspect that the discrepancy is due to the neglect of second 
(and higher) order corrections to the deviations from the guiding center motion in the analytic 
theory in Section 2. 

Superposed on the long-term periapse precessions are periodic variations on the same timescales 
(Ri 400 and 450 days for P2 and PI, respectively) as the epicyclic eccentricity variations seen in 
Fig. 3. For P2, the amplitude of the short-term variation is sufficiently large that the precession is 
retrograde when the eccentricity is large and prograde and faster than the long-term rate when the 
eccentricity is small. The short-term periodic variations in the epicyclic eccentricities and periapse 
longitudes are due to interactions between P2 and PI, and we can explain the periods by noting 
that their orbits are close to the 3:2 mean-motion commensurability (azimuthal period ratio = 1.539 
for the numerical integration, in agreement with the ratio 1.537 for the best-fit orbital periods in 
Table 1). The disturbing potential for the interactions between P2 and PI, which is not included 
in the analytic theory in Section 2, can be expanded into a cosine series in the usual manner (see, 
e.g., Murray and Dermott 1999), and because of the proximity to the 3:2 commensurability, we 
expect the interactions to be dominated by the terms in the disturbing potential associated with 
the cosine arguments (or resonance variables) 62 = 202 — 30i + VJ2 and 61 = 2(p2 — 30i + wi Fig. 
5 shows the evolution of 62 and 9i for 10^ days in the high albedo case. It is clear from Fig. 5 
that neither 62 nor 61 is in resonance and that the circulation periods of 62 and 9i are 400 and 
450 days, respectively. Thus the short-term periodic variations in 62 and W2 are associated with 
the circulation of 02 and those in ei and wi are associated with the circulation of 9i. 

The amplitudes of the short-term variations in e and m are larger for both P2 and PI for larger 
masses of P2 and PI, as long as the masses do not exceed about half the low-albedo ones. In the 
run with half the low-albedo masses (Figs. 6 and 7), ei and zui show significant variations, with the 
minimum ei much smaller than the initial ei = 0.0052. Both W2 and wi alternate between nearly 
linear retrograde precession when the eccentricities are large and very fast prograde precession when 
the eccentricities are small. However, the long-term trends of both zu2 and wi remain prograde, 
with periods of 1900 and 5100 days, respectively, which are slightly shorter than in the case with 
the high-albedo masses. 

In the run with the low-albedo masses (Figs. 8 and 9), while the properties of the orbit of P2 
follow the same trends with increasing masses as above (i.e., larger amplitudes for the short-term 
variations of 62 and ZU2 and a shorter period of 1800 days for the long-term prograde precession of 
W2), those of PI do not follow the previous trends. The amplitude of the variation in ei is now 
smaller than in the case with half the low-albedo masses, and wi shows a completely retrograde 



*For convenience, we define the resonance variables $2 and 9i using the true longitudes 02 and 0i instead of mean 
longitudes. While it is possible to define a mean longitude via a mean anomaly that is proportional to not, it cannot 
be easily computed from the position and velocity. 
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precession with a period of 500 days. The explanation can be found in the plot of the evolution of 

the resonance variables 62 and 9i in Fig. 10. While 62 circulates as in the high albedo case shown 
in Fig. 5 (the nearly empty region about 180° is due to the fast progradc precession of W2 when 62 
is small), 9i librates about 180°. Thus the case with the low-albedo masses differs from the cases 
with lower masses in having 61 in resonance. 

4. DISCUSSION AND CONCLUSIONS 

We have analyzed the orbits of the recently discovered satellites of Pluto, S/2005 PI and 
S/2005 P2. Because of the rather large mass ratio of Charon-Pluto {rric/mp ~ 0.1), the orbits of P2 
and PI are non-Keplerian even if P2 and PI have negligible masses. The analytic theory in Section 
2 with P2 and PI treated as test particles shows that the motion in R and (j) can be described by the 
superposition of the circular motion of the guiding center at mean motion riQ, the epicyclic motion 
represented by eccentricity e at epicyclic frequency kq, and the forced oscillations at frequencies 
k\nQ — Hpcl {k = 1, 2, 3, . . .) due to the non-axisymmetric components of the gravitational potential 
rotating at the mean motion Upc of Pluto-Charon; the motion in z at the vertical frequency uq 
decouples from that in R and (p. With i'q > uq > uk > kq, where uk is the Keplcrian mean motion 
about the center of mass of Pluto and Charon, we found that the azimuthal period Pq = 27r/no is 
shorter than the Keplcrian orbital period and that the periapse and ascending node (relative to the 
Pluto-Charon orbital plane) precess at nearly equal rates in opposite directions (prograde for the 
periapse and retrograde for the node). We have also performed a series of direct numerical orbit 
integrations with different assumed masses for P2 and PI, and the results presented in Section 3 
show the increasing effects of the proximity of the orbits of P2 and PI to the 3:2 mean-motion 
commensurability with increasing masses. As shown in Fig. 1, the deviation from Kepler's third 
law is already detected in the unperturbed Keplcrian fit of BGYYS (which was previously pointed 
out by BGYYS as discrepancies in the total mass of Pluto-Charon inferred from the orbits of 
Charon, P2, and PI). Since the other non-Keplerian behaviors depend on the masses of P2 and PI, 
a dynamical fit to the data that accounts for the interactions among Charon, P2, and PI should 
allow us to place constraints on the masses of P2 and PI, although the existing data consisting 
of only 12 observations over a 1-year time span may not be sufficient. (BGYYS included in their 
unperturbed Keplcrian fits the discovery data for P2 and PI from Weaver et al. 2006, which were 
taken two years after the last observation of BGYYS, but the discovery data have larger errors and 
may not be useful in constraining parameters other than the orbital periods.) 

If the albedos of P2 and PI are high and of order of that of Charon, the masses of P2 and PI 
are sufficiently low that their orbits are well described by the analytic theory. The largest correction 
due to the proximity to the 3:2 commensurability is a significant variation of the precession of ^2 

on the period of circulation (w 400 days) of the resonance variable O2 = 2<p2 — + W2 (Figs. 4 
and 5). However, the variation in the orbital radius R2 of P2 due to the forced oscillations are 
large (Fig. 2), and the eccentricity with a large error (0.0023 it 0.0021) found by the Keplcrian fit 
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of BGYYS probably results from an attempt to fit the 12 data points without taking into account 

the forced oscillations. Thus, there is at present no evidence that P2 has any significant epicyclic 
eccentricity, and it is likely that the periapsc precession of P2 would be difficult to measure. On 
the other hand, the orbit of PI has a significant epicyclic eccentricity. Over the 1-ycar time span 
of the existing data, the prograde precession of wi with a period of 5300 days (Fig. 4) would have 
resulted in a 25° change in wi, which may be difficult to detect with only 12 data points. But 
there would be more than 3.6 yr between the first existing data point and any additional data of 
comparable quality taken after the writing of this paper (Feb. 2006), and the > 90° precession of 
wi should be detectable. 

If the albedos of P2 and PI are low and of order of that of comets, the masses of P2 and PI 
are sufficiently large that there are significant short-term variations in their epicyclic eccentricities 
and/or periapse longitudes due to the proximity to the 3:2 commensurability. In the case with 
half the low-albedo masses (Figs. 6 and 7), there are significant variations in 62 and W2 on the 
circulation period of 400 days of 62, and in ei and tui on the circulation period of 450 days of 
9i = 2(j)2 — 3(l)i+wi. In the case with low-albedo masses (Figs. 8, 9, and 10), 9i is in resonance and 
zui shows a retrograde precession with a period of only 500 days. Since the existing data can be 
reasonably fitted by unperturbed Kcplcrian orbits, one might speculate that the existing data are 
already inconsistent with masses of P2 and PI near the upper end of the expected range. But we 
cannot rule out the possibility that orbits like those in Figs. 6-9 can be fitted by Keplerian orbits, 
if the orbits are sparsely sampled as in the existing data set. Another indication that masses of P2 
and PI near the upper end of the expected range may already be rnled out comes from the orbital 
eccentricity, Cc, of Charon. As pointed out by Stern et al. (1994), the perturbations from additional 
satellites in the Pluto-Charon system induce eccentricity in Charon's orbit, and the observed value 
or upper limit of Cc can be used to constrain the masses of the additional satellites. In Fig. 11 we 
show the variation of Cc for 800 days in the run with low-albedo masses. The eccentricity Cc varies 
up to 2 X 10~^, which is significantly larger than the best-fit Cc = 0.0 ± 7.0 x 10~^ in Table 1. But 
it should be noted that the random error of 7.0 x 10~^ corresponds to shifts of the order of 0.1 
mas in the positions of Charon relative to Pluto, and it is unclear that systematic error in, e.g., the 
correction from the center of light to the center of body is below 0.1 mas. In any case, a data set 
that samples more densely any possible variations in the orbits of P2 and PI on the 400-500-day 
timescales should provide strong constraints on the masses of P2 and PI. 

There are no apparent effects due to the respective proximity of P2 and PI to the high-order 
4:1 and 6:1 mean-motion commensurabilities with Charon. However, Ward and Canup (2006) have 
proposed that P2 and PI may have been trapped in the corotation resonances at these commensu- 
rabilities (i.e., those associated with the resonance variable (pc — 4:(j)2+3^c for P2 and (/>c — 6(/)i -|- SzUc 
for PI) during the tidal expansion of Charon's orbit, if Charon formed with large orbital eccentricity 
from a giant impact on Pluto. P2 and PI escaped from the corotation resonances when Charon's 
orbital eccentricity was tidally damped to very small value. 

Our analysis shows that continued observations of the Pluto system with HST (and possibly 
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ground-based adaptive optics) in the near future will allow us to detect the non-Keplerian behaviors 
of the orbits of P2 and PI (in addition to the already detected deviation from Kepler's third law) 
and to thereby constrain their masses. Much more precise determination of the orbits and masses 
of P2 and PI will be possible as the New Horizons spacecraft approaches the Pluto system in 2015. 

We thank Robin Canup and Alan Stern for sending us preprints on the new satellites of Pluto. 
This research was supported in part by NASA grant NNG05GK58G. 
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Table 1. Orbital Parameters at Epoch JD 2452600.5 from Keplerian fits by Buie et al. (2006) 



Parameter 


Charon 


S/2005 P2 


S/2005 PI 


Period P (days) 


6.3872304(11) 


24.8562(13) 


38.2065(14) 


Semimajor Axis a (km) 


19571.4(4.0) 


48675(121) 


64780(88) 


Eccentricity e 


0.000000(70) 


0.0023(21) 


0.0052(11) 


Inclination i (deg) 


96.145(14) 


96.18(22) 


96.36(12) 


Long. Ascending Node 17 (deg) 


223.046(14) 


223.14(23) 


223.173(86) 


Long. Periapse vj (deg) 




216(13) 


200.1(3.7) 


Mean Long, at Epocfi A (deg) 


257.946(13) 


123.14(20) 


322.71(23) 



Note. — The parameters are for the orbit of Charon relative to Pluto and the 
orbits of PI and P2 relative to the center of mass of Pluto-Charon. Numbers 
in parentheses are 1 a errors in the least significant digits. 
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Table 2. Parameters of Analytic Theory 



Parameter 


S/2005 P2 


S/2005 PI 


Ro (km) 


48675 


64780 


Pk = 27r/nK (days) 


25.0518 


38.4628 


no/ UK 


1.00635 


1.00341 




0.99198 


0.99612 


vq/tik 


1.02053 


1.01063 


Ci 


0.001275 


0.000149 


C2 


0.001373 


0.000228 


Cs 


0.000204 


0.000026 


C4 


0.000044 


0.000004 


Di 


0.003220 


0.000458 


D2 


0.006813 


0.001764 


D3 


0.001496 


0.000314 


L>4 


0.000437 


0.000072 
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Fig. 1. — Deviation of azimuthal period Pq from Keplerian orbital period Pk as a function of the 
guiding center radius Rq. The solid line shows Pr/Pq — 1 = no/riK — 1 from Eq. (17) for the best- 
fit mass ratio of Charon-Pluto, mc/mp = 0.1165, while the dotted lines show the same quantity 
for m^c/nip that are la (±0.0055) from the best-fit value. The values of Pk/P — 1 and a for the 
satellites P2 and PI, with orbital periods P and semimajor axes a from unperturbed Keplerian fits 
(Table 1), are shown with their la error bars. 
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Fig. 2. — Variations in the orbital radii R2 and Ri of P2 and PI, respectively, for 800 days in the 
numerical orbit integration with the high-albcdo masses for P2 and PI. The dashed lines indicate 
the semimajor axes, a, and the maximum and minimum radii, a(l it e), for the best-fit Keplerian 
orbits in Table 1. The variation in R2 is dominated by high-frequency forced oscillations due to the 
non-axisymmetric components of the potential rotating at the mean motion of Pluto-Charon, but 
a small epicyclic motion with period 25.2 days is also visible. The variation in Ri is dominated 
by the epicyclic motion with eccentricity 0.0052 and period 38.6 days. 
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Fig. 3. — Same as Fig. 2, but for the transformed orbital radii i?2 and R[ (Eq. [35]). Most of 
the forced oscillations in R2 and Ri are eliminated, and small periodic variations in the maximum 
and minimum values of i?2 and R[ , or equivalently in the epicyclic eccentricities 62 and ei , become 
visible. The periods are about 400 days for 62 and 450 days for ei. 
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Fig. 4. — Evolution of the longitudes of pcriapse W2 and zui of P2 and PI, respectively, for 
10^ days in the high albedo case. There are short-term periodic variations on the same timescales 
as 62 and ei , but the long-term evolution is prograde precessions with periods of 2000 days for P2 
and 5300 days for PI. 
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Fig. 5. — Same as Fig. 4, but for the resonance variables 62 = 2^2— 30i+tU2 and 9i = 2^2— 3^i+tui 
at the 3:2 commensurabihty between P2 and PI. The period of the short-term variations in 62 and 
W2 (ei and wi) is the circulation period of 62 (Gi)- 
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Fig. 6. — Same as Fig. 2, but for the numerical integration with P2 and PI having half the low- 

albcdo masses. The periodic variations in 62 and ei have amplitudes that are significantly larger 
than in the high-albedo case (Fig. 3) and are visible without transforming R2 and Ri to R'2 and 
R',. 
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Fig. 7. — Same as Fig. 4, but for the numerical integration with P2 and PI having half the 
low-albedo masses. Compared to the high albedo case (Fig. 4), the amplitudes of the short-term 
periodic variations are significantly larger, and the periods of the long-term prograde precession 
are slightly shorter. 
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Fig. 8. — Same as Fig. 2, but for the numerical integration with the low-albedo masses for P2 and 
PI. The amplitude of the periodic variation in 62 continues to increase with mass, but that in ei 
is smaller than in the case with half the low-albedo masses (Fig. 6). 
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Fig. 9. — Same as Fig. 4, but for the numerical integration with the low-albcdo masses for P2 
and PI. While vd2 continues to show long-term prograde precession with large short-term periodic 
variation, w\ shows a completely retrograde precession with a period of 500 days. 
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Fig. 10. — Same as Fig. 5, but for the numerical integration with the low-albcdo masses for P2 
and PI. While Q2 circulates as in the high albedo case shown in Fig. 5 (the nearly empty region 
about 180° is due to the fast prograde precession of W2 when 62 is small), Q\ librates about 180°. 
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Fig. 11. — Variation in the orbital eccentricity Cc of Charon for 800 days in the numerical integration 
with the low-albedo masses for P2 and PI. The dashed line indicates Cc = 7.0 x 10~^, which is 1 a 
above the best-fit value (ec = 0.0) in Table 1. 



